//==============================================================================
//
//   DRndDist - the distributed pseudo random generator in the ofc-library
//
//               Copyright (C) 2007  Dick van Oudheusden
//  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public
// License along with this library; if not, write to the Free
// Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//
//==============================================================================
// 
//  $Date: 2008-08-03 18:31:55 $ $Revision: 1.6 $
//
//==============================================================================

/*
 * This class is based on the algorithms in the gsl library.
 *
 */

#include <math.h>

#include "ofc/DRndDist.h"


#if _INTERFACE_

#include <stdint.h>

#include "ofc/config.h"

#include "ofc/DRandomable.h"


//
// The DRndDist class implements pseudo random generators distributions.
//
// @example
// #include <stdio.h>
// #include "ofc/DRndDist.h"
// #include "ofc/DMRnd.h"
// 
// int main(int argc, char *argv[])
// {
//   DMRnd *rnd = [DMRnd new];
//   
//                                        // Uniform distributed random number
//   printf("Random number, uniform distributed, range 0.0..10.0    : %f.\n", 
//     [DRndDist nextUniform :rnd :0.0 :10.0]);
// 
//                                        // Normal distributed random number
//   printf("Random number, normal distributed,  range 0.0..10.0    : %f.\n",
//     [DRndDist nextNormal :rnd :0.0 :10.0]);
// 
//                                        // Exponential distributed random number
//   printf("Random number, exponential distributed, mean 5.0       : %f.\n",
//     [DRndDist nextExponential :rnd :5.0]);
// 
//                                        // Gamma distributed random number
//   printf("Random number, gamma distributed, alpha 2.0, beta 3.0  : %f.\n",
//     [DRndDist nextGamma :rnd :2.0 :3.0]);
// 
//                                        // Beta distributed random number
//   printf("Random number, beta distributed, alpha 2.0, beta 3.0   : %f.\n",
//     [DRndDist nextBeta :rnd :2.0 :3.0]);
// 
//                                        // Binomial distributed random number
//   printf("Random number, binomial distributed, prob .5 runs 10   : %u.\n",
//     [DRndDist nextBinomial :rnd :0.5 :10]);
// 
//                                        // Poisson distributed random number
//   printf("Random number, poisson distributed, mu 5.0             : %u.\n",
//     [DRndDist nextPoisson :rnd :5.0]);
// 
//                                        // Weibull distributed random number
//   printf("Random number, weibull distributed, alpha 2.0, beta 3.0: %f.\n",
//     [DRndDist nextWeibull :rnd :2.0 :3.0]);
// 
//   [rnd free];                          // Cleanup
// 
//   return 0;
// }
//

@interface DRndDist : Object
{
}

#endif


@implementation DRndDist

  
//// Random Generator Distribution class methods

//
// Generate a positive uniform random number <0,1>
// 
// @param rnd      the uniform random genetor
// 
// @param the positive uniform random number
// 
+ (double) nextPosDouble :(id <DRandomable>) rnd
{
  double r;
  
  do
  {
    r = [rnd nextDouble];
  }
  while (r == 0.0);
  
  return r;
}


//
// Generate a random number with a uniform distribution in the range of [start,end>
// 
// @param rnd      the uniform random generator
// @param start    the start of the range
// @param end      the end of the range
//
// @return the random number
//

+ (double) nextUniform :(id <DRandomable>) rnd  :(double) start :(double) end
{
  return [rnd nextDouble :start :end];
}


//
// Generate a random number with a normal or gaussian distribution
// 
// @param rnd      the uniform random generator
// @param mu       the mean
// @param sigma    the standard deviation
//
// @return the random number
//

+ (double) nextNormal :(id <DRandomable>) rnd  :(double) mu :(double) sigma
{
  /* Ratio method (Kinderman-Monahan); see Knuth v2, 3rd ed, p130.
   * K+M, ACM Trans Math Software 3 (1977) 257-260.
   * With Leva's modifications to the original K+M method; see:
   * J. L. Leva, ACM Trans Math Software 18 (1992) 449-453 and 454-455. */
  
  double u, v, x, y, Q;

  do
  {
    u = 1 - [rnd nextDouble];
    v = [rnd nextDouble] - 0.5;

    v *= 1.7156;

    x = u - 0.449871;
    y = fabs (v) + 0.386595;
    Q = x * x + y * (0.19600 * y - 0.25472 * x);
  }
  while ((Q >= 0.27597) && ((Q > 0.27846) || (v * v > -4 * u * u * log (u))));

  return mu + sigma * (v / u);
}


//
// Generate a random number with an exponential distribution
// 
// @param rnd      the random generator
// @param mu       the mean (>0)
//
// @return the random number
//

+ (double) nextExponential :(id <DRandomable>) rnd :(double) mu
{
  if (mu <= 0.0)
  {
    WARNING(DW_INVALID_ARG, "mu");
    
    return 0.0;
  }
  else
  {
    double u = [DRndDist nextPosDouble :rnd];

    return -mu * log (u);
  }
}


//
// Generate a random number with a gamma distribution
// 
// @param rnd      the uniform random generator
// @param alpha    the alpha parameter (>0)
// @param beta     the beta parameter (>0)
//
// @return the random number
//

+ (double) nextGamma :(id <DRandomable>) rnd :(double) alpha :(double) beta
{
  if ((alpha <= 0.0) || (beta <= 0.0))
  {
    WARNING(DW_INVALID_ARG, "alpha/beta");
    
    return 0.0;
  }
  else
  {
    /* Based on Marsaglia and Tsang, "A Simple Method for
     * generating gamma variables", ACM Transactions on Mathematical
     * Software, Vol 26, No 3 (2000), p363-372.
     */
    
    if (alpha < 1)
    {
      double g = [self nextGamma :rnd :1.0 + alpha :beta];
      double u = [DRndDist nextPosDouble :rnd];
    
      return g * pow (u, 1.0 / alpha);
    }
    else
    {
      double x, v, u;
      double d = alpha - 1.0 / 3.0;
      double c = (1.0 / 3.0) / sqrt (d);

      while (1)
      {
        do
        {
          x = [DRndDist nextNormal :rnd :0.0 :1.0];
          v = 1.0 + c * x;
        }
        while (v <= 0);

        v = v * v * v;
        u = [DRndDist nextPosDouble :rnd];

        if (u < 1 - 0.0331 * x * x * x * x) 
          break;

        if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
          break;
      }
    
      return beta * d * v;
    }
  }
}


//
// Generate a random number with a beta distribution
// 
// @param rnd      the uniform random generator
// @param alpha    the alpha parameter (>0)
// @param beta     the beta parameter (>0)
//
// @return the random number
//

+ (double) nextBeta :(id <DRandomable>) rnd :(double) alpha :(double) beta
{
  if ((alpha <= 0.0) || (beta <= 0.0))
  {
    WARNING(DW_INVALID_ARG, "alpha\beta");
    
    return 0.0;
  }
  else
  {
    /* Based on Knuth */

    double x1 = [self nextGamma :rnd :alpha :1.0];
    double x2 = [self nextGamma :rnd :beta  :1.0];

    return x1 / (x1 + x2);
  }
}


//
// Generate a random number with a binomial distribution
// 
// @param rnd      the uniform random generator
// @param p        the probability [0,1]
// @param n        the number of trails (0>)
//
// @return the random number
//

+ (unsigned int) nextBinomial :(id <DRandomable>) rnd :(double) p :(unsigned int) n
{
  unsigned int k = 0;

  if ((p < 0.0) || (p > 1.0))
  {
    WARNING(DW_INVALID_ARG, "p");
  }
  else
  {
    /* Based on Knuth */

    unsigned int i, a, b = 0;

    while (n > 10)
    {
      double X;
    
      a = 1 + (n / 2);
      b = 1 + n - a;

      X = [DRndDist nextBeta :rnd :(double) a :(double) b];

      if (X >= p)
      {
        n = a - 1;
        p /= X;
      }
      else
      {
        k += a;
        n = b - 1;
        p = (p - X) / (1 - X);
      }
    }

    for (i = 0; i < n; i++)
    {
      double u = [rnd nextDouble];
    
      if (u < p)
      {
        k++;
      }
    }
  }
  
  return k;
}


//
// Generate a random number with a poisson distribution
//  
// @param rnd      the uniform random generator
// @param mu       the mean (>=0)
//
// @return the random number
//

+ (unsigned int) nextPoisson :(id <DRandomable>) rnd :(double) mu
{
  if (mu < 0.0)
  {
    WARNING(DW_INVALID_ARG, "mu");
    
    return 0.0;
  }
  else
  {
    /* Based on Knuth */

    double   emu;
    double   prod = 1.0;
    unsigned 
    int      k    = 0;

    while (mu > 10)
    {
      unsigned int m = mu * (7.0 / 8.0);

      double X = [DRndDist nextGamma :rnd :(double) m :1.0];

      if (X >= mu)
      {
        return k + [DRndDist nextBinomial :rnd :(mu / X) :(m - 1)];
      }
      else
      {
        k  += m;
        mu -= X; 
      }
    }


    emu = exp (-mu);

    do
    {
      prod *= [rnd nextDouble];
      k++;
    }
    while (prod > emu);

    return k - 1;
  }
}


//
// Generate a random number with a pareto distribution
//  
// @param rnd      the uniform random generator
// @param alpha    the scale parameter (>0)
// @param beta     the shape parameter (>0)
//
// @return the random number
//

+ (double) nextPareto :(id <DRandomable>) rnd :(double) alpha : (double) beta
{
  if ((alpha <= 0.0) || (beta <= 0.0))
  {
    WARNING(DW_INVALID_ARG, "alpha\beta");
    
    return 0.0;
  }
  else
  {
    double x = [DRndDist nextPosDouble :rnd];

    double z = pow (x, -1 / beta);

    return alpha * z;
  }
}


//
// Generate a random number with a weibull distribution
//  
// @param rnd      the uniform random generator
// @param alpha    the scale parameter (0>)
// @param beta     the shape parameter (0>)
//
// @return the random number
//

+ (double) nextWeibull :(id <DRandomable>) rnd :(double) alpha :(double) beta
{
  if ((alpha <= 0.0) || (beta <= 0.0))
  {
    WARNING(DW_INVALID_ARG, "alpha\beta");
    
    return 0.0;
  }
  else
  {
    double x = [DRndDist nextPosDouble :rnd];

    double z = pow (-log (x), 1 / beta);

    return alpha * z;
  }
}

@end

/*==========================================================================*/
